3D DEM Simulations and Experiments on Spherical Impactor Penetrating into the Elongated Particles

In this study, a brass or glass spherical impactor vertically penetrating into a granular bed composed of mono-sized spherical or elongated particles was simulated with three-dimensional (3D) discrete element method (DEM). Good agreement of the particle masses in the cup before and after penetration can be found in the simulations and experiments. The effects of particle length (Lp), friction coefficient, and particle configuration on the penetration depth of the impactor, ejecta mass, and solid volume fraction describing the response of the granular bed are discussed. The penetration depth is negatively correlated with Lp as the corresponding solid volume fraction of the granular bed decreases. A smaller friction coefficient leads to a larger penetration depth of the impactor and more ejection of particles. When the impactor is penetrating the Lp = 10 mm elongated particles, the penetration depth is negatively correlated to the order parameter and solid volume fraction.


Introduction
Non-spherical particles are widely ubiquitous in nature and industry; for instance, micro-particles of various shapes have been applied in drug transportation and mixing. Particle shape has been shown to influence various processes including drug transport and crystal preparation [1][2][3][4]. When non-spherical elongated particles are packed, the internal microstructure of the granular bed has a complex and important influence on the macroscopic mechanical response of the material, which is different from that of spherical particles. Furthermore, the mechanical behavior of non-spherical particles affected by the impact will have a significant influence on transportation and industrial production. Therefore, the study of the impact into elongated particles is of great importance.
The finite element method (FEM) is an effective numerical method for continuum material, such as variational phase-field problems [5,6], quasistatic frictional contact problems [7], and elastic systems [8]. In addition, the discrete element method (DEM) is one of the most effective numerical methods for studying the dynamic properties of granular medium; it is a discontinuous numerical method first proposed by Cundall in 1971 [9]. Based on the theory of Newton's second law and the contact model between particles, the DEM can be applied to the particle dynamics and kinematics analysis by calculating the displacement and force of each particle in the granular system based on the explicit time step iteration.
Both 2D [10][11][12][13][14] and 3D [15][16][17] discrete element models have been used in numerical simulations of object impacts into target granular medium. Studies have shown that the response of the granular bed and impactor depends on many factors, such as the shape, size, angle, and penetration velocity of the impactor [10,[18][19][20][21], and the friction coefficient and porosity of the granular bed [11]. In addition, the energy dissipates drastically during the transient impact process, which has a complex relationship with the physical parameters or material properties of the granular medium [22]. Based on studies of spherical particles, the effect of elongated particles can be discussed comparatively.
Numerical models of regular non-spherical particles such as cylinders [23], spherocylinders [24,25], polyhedrons [26,27], super-ellipsoids [28,29], multi-super-ellipsoids [30,31], and arbitrarily shaped elements [32] have been developed gradually. The packing of elongated particles was discussed by changing various parameters, such as coefficient of friction and coefficient of restitution [33]. The behavior of elongated particles with varying lengths was explored by discharging them from a rectangular hopper [34]. Studies have shown that the length of the particle can be used to adjust the buffer capacity for the non-spherical particle system [35]. Besides, the effects of particle length and configuration of granular bed on impact have not been discussed, and the effect mechanisms of the elongated particle system have not been studied comprehensively. This study can provide a fundamental insight for the mechanical behavior of impact into non-spherical particles.
In this study, we carried out experiments on vertical impact into granular material and compared the ejecta masses and initial solid volume fractions of the granular bed between experiments and simulations, considering spherical and elongated particles as the objective granular medium. For the simulations based on DEM, we investigated the behaviors of the ejected particles, the particles in the cup, and the impactor. For the above three objects of study, we considered the effects of some key factors according to their different behavior characteristics, including the length of elongated particles, the friction coefficient, the shape of the granular cup, and the configuration of particles.

Numerical Model and Methodology
In this study, the translational and rotational movements of particles are governed by Newton's second law of motion. The equations are written as follows: in which m i , v i , I i , and ω i are the mass, translational velocity, moment of inertia, and rotational velocity of particle i, respectively, and F i and T i are external forces and torques exerting on it. The contact types of the sphero-cylindrical particle for the DEM simulation are classified into four groups according to the normal contact forces, as shown in Table 1. The normal contact force calculation and implemented contact detection in this study were proposed by Kidokoro et al. [25] and modified by Guo et al. [36]. The magnitudes and directions of normal and tangential contact forces are determined by contact position and the overlap of the two contacting particles. The tangential force model is the Mindlin model, which only takes into account the static process in this study.

Scenarios Models
Group Ⅰ F n = 4 in which δn is the overlap in normal direction, E * is equivalent Young's modulus, and R * is equivalent radius of two objects in contact, which are defined as in which α is determined by the shape of contact area and determined to be 0.974 here.

Scenarios Models
Group II Table 1. Contact force models used in the DEM simulations.

Scenarios Models
Group Ⅰ F n = 4 3 * * (3) in which δn is the overlap in normal direction, E * is equivalent Young's modulus, and R * is equivalent radius of two objects in contact, which are defined as in which α is determined by the shape of contact area and determined to be 0.974 here.
Group Ⅲ = ⎩ ⎨ in which θ is the angle between the two major axes of particles. = √ √ * * , when θ is equal to zero.
is the same as Equation (6).
in which k is a constant determined as 2.5, l is the length of contact area along the major axis, b is the width of contact area, = 2 * .
Tangential force model where Ft 0 and Ft are the tangential force vectors in the previous time step and the current time step, respectively. G * is governed by * = + , in which G1 and G2 are the shear moduli of the two objects in contact, ν1 and ν2 are the corresponding Poisson's ratios. is the effective radius of contact, = * , and dt represents the incremental tangential displacement in the present time step.
Note: The equations are summarized according to the proposed models [36].
For the normal and tangential damping force [36], F n d and F t d are determined by where c is equal to 1 if the normal contact force Fn is proportional to , and c is equal to if the normal contact force Fn is proportional to . m * is equivalent mass, m * = in which α is determined by the shape of contact area and determined to be 0.974 here.

Scenarios Models
Group Ⅰ F n = 4 in which δn is the overlap in normal direction, E * is equivalent Young's modulus, and R * is equivalent radius of two objects in contact, which are defined as in which α is determined by the shape of contact area and determined to be 0.974 here.
Group Ⅲ = ⎩ ⎨ in which θ is the angle between the two major axes of particles. = √ √ * * , when θ is equal to zero.
is the same as Equation (6).
in which k is a constant determined as 2.5, l is the length of contact area along the major axis, b is the width of contact area, = 2 * .
Tangential force model where Ft 0 and Ft are the tangential force vectors in the previous time step and the current time step, respectively. G * is governed by * = + , in which G1 and G2 are the shear moduli of the two objects in contact, ν1 and ν2 are the corresponding Poisson's ratios. is the effective radius of contact, = * , and dt represents the incremental tangential displacement in the present time step.
Note: The equations are summarized according to the proposed models [36].
For the normal and tangential damping force [36], F n d and F t d are determined by where c is equal to 1 if the normal contact force Fn is proportional to , and c is equal to if the normal contact force Fn is proportional to . m * is equivalent mass, m * = in which θ is the angle between the two major axes of particles.
n , when θ is equal to zero. F max n is the same as Equation (6).

Scenarios Models
Group Ⅰ F n = 4 in which δn is the overlap in normal direction, E * is equivalent Young's modulus, and R * is equivalent radius of two objects in contact, which are defined as in which α is determined by the shape of contact area and determined to be 0.974 here.
Group Ⅲ = ⎩ ⎨ in which θ is the angle between the two major axes of particles. = √ √ * * , when θ is equal to zero.
is the same as Equation (6).
in which k is a constant determined as 2.5, l is the length of contact area along the major axis, b is the width of contact area, = 2 * .
Tangential force model where Ft 0 and Ft are the tangential force vectors in the previous time step and the current time step, respectively. G * is governed by * = + , in which G1 and G2 are the shear moduli of the two objects in contact, ν1 and ν2 are the corresponding Poisson's ratios. is the effective radius of contact, = * , and dt represents the incremental tangential displacement in the present time step.
Note: The equations are summarized according to the proposed models [36].
For the normal and tangential damping force [36], F n d and F t d are determined by where c is equal to 1 if the normal contact force Fn is proportional to , and c is equal to if the normal contact force Fn is proportional to . m * is equivalent mass, m * =  (6) in which k is a constant determined as 2.5, l is the length of contact area along the major axis, b is the width of contact area, b = √ 2R * δ n .
Tangential force model where F t 0 and F t are the tangential force vectors in the previous time step and the current time step, respectively. G * is governed by 1 G * = 2−ν 1 G 1 + 2−ν 2 G 2 , in which G 1 and G 2 are the shear moduli of the two objects in contact, ν 1 and ν 2 are the corresponding Poisson's ratios. a is the effective radius of contact, a = √ R * δ n , and v t c dt represents the incremental tangential displacement in the present time step.
Note: The equations are summarized according to the proposed models [36].
For the normal and tangential damping force [36], F d n and F d t are determined by where c is equal to 1 if the normal contact force F n is proportional to δ n , and c is equal to 5 6 if the normal contact force F n is proportional to δ n 3 2 . m * is equivalent mass, , and m 1 and m 2 are the masses of two contacting particles. v n and v t are the normal and tangential components of relative velocity, respectively. β = −lne √ π 2 +(lne) 2 is the contact damping coefficient, where e is the coefficient of restitution. S n is normal contact stiffness, and S t is tangential contact stiffness, given by S n = dF n dδ n and S t = 8G * a, respectively.

Experimental and Numerical Setup
The experimental apparatus used in this study is the same as what we used in the previous study [17]. Figure 1 shows the schematic of the experimental apparatus. The experiment was carried out in an airtight vacuum chamber (<200 mTorr), taking no account of the gas drag on the particles during the impact process. The impactor has a steel dot pasted on the top. Firstly, the impactor was fixed by an electromagnet at the central position of the crossbar at a height of H = 661 mm above the top of the granular bed. Then, the cylindrical particles were poured into the cup, and the surface of the granular bed was smoothed with a scraper. The total mass of the cup and the particles in it, m 0 , was weighed. The air was then pumped out until the pressure in the chamber was lower than 200 mTorr. Finally, the impactor was released from static state under gravity, and the vertical velocity of the impactor was calculated by V 0 = 2gH (where g is the acceleration of gravity). After all particles reached the stable state, the total mass of the cup and particles remaining in it, m, was weighed. The mass of ejected particles can be obtained by the equation ∆m = m 0 − m. The objective particles used in the experiments are the steel balls and steel cords cut into different lengths (L p = 4 mm and 6 mm), whose diameters are 2 mm. The impactors are brass and glass spheres, whose impact velocities are 3.60 m/s. The parameters of all the materials used in experiment are shown in Table 2.
pasted on the top. Firstly, the impactor was fixed by an electromagnet at the central position of the crossbar at a height of H = 661 mm above the top of the granular bed. Then, the cylindrical particles were poured into the cup, and the surface of the granular bed was smoothed with a scraper. The total mass of the cup and the particles in it, m0, was weighed. The air was then pumped out until the pressure in the chamber was lower than 200 mTorr. Finally, the impactor was released from static state under gravity, and the vertical velocity of the impactor was calculated by = 2g (where g is the acceleration of gravity). After all particles reached the stable state, the total mass of the cup and particles remaining in it, m, was weighed. The mass of ejected particles can be obtained by the equation Δm = m0 − m. The objective particles used in the experiments are the steel balls and steel cords cut into different lengths (Lp = 4 mm and 6 mm), whose diameters are 2 mm. The impactors are brass and glass spheres, whose impact velocities are 3.60 m/s. The parameters of all the materials used in experiment are shown in Table 2. The setup of DEM simulations is shown in Figure 2. The geometric and physical parameters of the cup, granular particles, and impactors are chosen from the experimental study. We mainly focus on the particles in the cup and the initial trajectory of ejected particles, and the cuboidal active simulation domain is big enough and will not affect the motion of ejected particles. The implementation of DEM simulation in this study was as follows: Firstly, the computing domain was divided into same-sized cuboid cells, and a certain number of particles were generated without contact and mapped into the cells. The particles were generated within the cylindrical domain with the same diameter as the cup and enough height. The generated particles then fell downward under gravity and packed in the cup until the cup was filled up with the particles. The translational and   The setup of DEM simulations is shown in Figure 2. The geometric and physical parameters of the cup, granular particles, and impactors are chosen from the experimental study. We mainly focus on the particles in the cup and the initial trajectory of ejected particles, and the cuboidal active simulation domain is big enough and will not affect the motion of ejected particles. The implementation of DEM simulation in this study was as follows: Firstly, the computing domain was divided into same-sized cuboid cells, and a certain number of particles were generated without contact and mapped into the cells. The particles were generated within the cylindrical domain with the same diameter as the cup and enough height. The generated particles then fell downward under gravity and packed in the cup until the cup was filled up with the particles. The translational and rotational motions of particles were controlled by Newton's second law of motion. The normal and tangential contact forces were calculated by the overlap between two contacting objects according to various contact models. When all particles settled down, the impactor was then released with an initial vertical velocity of 3.60 m/s right above the particle bed and penetrated into the particle bed in the center. rotational motions of particles were controlled by Newton's second law of motion. The normal and tangential contact forces were calculated by the overlap between two contacting objects according to various contact models. When all particles settled down, the impactor was then released with an initial vertical velocity of 3.60 m/s right above the particle bed and penetrated into the particle bed in the center. In the DEM simulations, the particles are steel spheres and elongated steel particles of Lp = 4 mm, 6 mm, 8 mm, 10 mm with a hemisphere cap at each end. The diameters of the spheres and elongated particles are equal to 2 mm, which is the same as the cylindrical particles used in experiments. Brass and glass spherical impactors with the same diameter are used in both experiments and simulations. The diameter, the height, the density, the coefficient of friction, and the coefficient of restitution are measured experimentally. Young's modulus and Poisson's ratio are selected from a handbook [37]. The physical and geometric parameters of particles and impactors are summarized in Table 2. Contact detection of the elongated particles in this study is the same as the sphero-cylinder model in the literature [24]. The difference is that the elongated sphero-cylinder used here is rigid, and it was a flexible elastic model in the literature.   In the DEM simulations, the particles are steel spheres and elongated steel particles of L p = 4 mm, 6 mm, 8 mm, 10 mm with a hemisphere cap at each end. The diameters of the spheres and elongated particles are equal to 2 mm, which is the same as the cylindrical particles used in experiments. Brass and glass spherical impactors with the same diameter are used in both experiments and simulations. The diameter, the height, the density, the coefficient of friction, and the coefficient of restitution are measured experimentally. Young's modulus and Poisson's ratio are selected from a handbook [37]. The physical and geometric parameters of particles and impactors are summarized in Table 2. Contact detection of the elongated particles in this study is the same as the sphero-cylinder model in the literature [24]. The difference is that the elongated sphero-cylinder used here is rigid, and it was a flexible elastic model in the literature.
The time step and boundary effect are two important parameters for DEM simulations; therefore, the sensitivities of these two parameters are studied before we determine the final values of them. The penetration depths of five different time steps are compared. The L p = 10 mm steel elongated particles and the glass impactor are chosen for simulation, and the parameters of materials are shown in Table 2.    Table 2).

Results and Discussion
In this section, the key parameters on impact are discussed, including particle length, friction coefficient, cup shape, and particle configuration. In addition, three objects will be studied, including the particles out of the cup, the particles in the cup, and the impactor. Their fundamental parameters are analyzed, such as the ejecta mass, the penetration depth, the solid volume fraction, and the kinetic energy.

Effect of Particle Length
The effect of particle length on impact is studied in this section, and spherical particles and elongated particles of four different lengths with the same diameter are discussed. Two materials of brass and glass impactor with impact velocities of 3.60 m/s are used. The friction coefficients are the same as the base case and material parameters are shown in Table 2.
The kinetic energy of the ejected particles, the impactor, and the particles left in the cup are calculated. The kinetic energy includes rotational kinetic energy and translational kinetic energy. For particle i with mass mi and moment of inertia Ii, the translational  Table 2).

Results and Discussion
In this section, the key parameters on impact are discussed, including particle length, friction coefficient, cup shape, and particle configuration. In addition, three objects will be studied, including the particles out of the cup, the particles in the cup, and the impactor. Their fundamental parameters are analyzed, such as the ejecta mass, the penetration depth, the solid volume fraction, and the kinetic energy.

Effect of Particle Length
The effect of particle length on impact is studied in this section, and spherical particles and elongated particles of four different lengths with the same diameter are discussed. Two materials of brass and glass impactor with impact velocities of 3.60 m/s are used. The friction coefficients are the same as the base case and material parameters are shown in Table 2.
The kinetic energy of the ejected particles, the impactor, and the particles left in the cup are calculated. The kinetic energy includes rotational kinetic energy and translational kinetic energy. For particle i with mass m i and moment of inertia I i , the translational velocity and rotational velocity at time t are defined as v i (t) and w i (t), respectively, and the total kinetic energy of N particles is given by The initial kinetic energy of the impactor is E k0 = 1 2 MV 0 2 , and the initial potential energy of the impactor is E p0 = MgH 0 (the position of zero gravitational potential energy is at the final stopping point of the impactor). The mass of impactor M is calculated by the density and size shown in Table 2. We take the percentage of the energy to the sum of E k0 and E p0 as the scaled energy. ticles is close to that of 4 mm, and their values are only 0.5%, indicating that the total mass of ejected particles is very small. Similarly, when the length of particles increases to 8 mm, the ejecta mass for the brass impactor is almost the same as that of 6 mm. By comparing the results of the two impactors, it can be concluded that when the mass of the elongated particle increases to a certain value, the transferred energy cannot support the ejection of more particles.   Table 2). Figure 4a is the dynamic change of ejecta mass and Figure 4b is the comparison of the final ejecta masses between simulations and experiments. In Figure 4a, we set the time that the impactor touches the top of the particle bed as t = 0 ms. The end time of particle ejection for brass and glass impactors is about 90 ms. It is obvious that the ejecta masses of the brass impactor are larger than those of the glass impactor due to higher impactor energy. For the brass impactor, the ejecta masses of L p = 6 mm, 8 mm, and 10 mm elongated particles are closer and smaller than that of L p = 4 mm. For the glass impactor, the ejecta masses of L p = 4 mm, 6 mm, 8 mm, and 10 mm elongated particles are approximately equal and smaller than that of 2 mm spheres.
In Figure 4b, the agreement of solid volume fractions and ejecta masses between experiments and simulations can be found. The results show that the ejecta mass of spherical particles is much larger than that of elongated particles. In this simulation, the method of changing the aspect ratio is increasing the length of particles (keeping the diameter unchanged). In this way, the mass of a single particle also increases as a result of increasing particle length. When the length of particles increases to 6 mm, the impact energy of the glass impactor cannot make more elongated particles eject from the cup because of the large mass of a single particle. Therefore, the percentage of total ejecta mass of 6 mm particles is close to that of 4 mm, and their values are only 0.5%, indicating that the total mass of ejected particles is very small. Similarly, when the length of particles increases to 8 mm, the ejecta mass for the brass impactor is almost the same as that of 6 mm. By comparing the results of the two impactors, it can be concluded that when the mass of the elongated particle increases to a certain value, the transferred energy cannot support the ejection of more particles. Figure 5a,b show the kinetic energy of ejected particles for two impactors. The E e k reached its peak value at about 2 ms when the particles collided with the impactor and gained the maximum velocities. Some particles began to launch upwards and fly in the air along parabolic trajectories. The ejected particles were moving upward along the parabolic trajectory before 20 ms, and their velocities were hence decreasing, resulting in the decreasing kinetic energy. It is easy to find that the peak values of E e k s for spherical particles are much larger than those of elongated particles for both glass and brass impactors because spherical particles roll more easily than elongated ones. Comparatively, the E e k s of different elongated particles lengths are obviously smaller than those of spherical particles, and E e k s are slightly smaller as the particles become longer, whether for the glass or brass impactor. We can suppose that the kinetic energy transferred to the ejected particles is to support their movement and more energy will be needed to drive a particle moving if its mass becomes larger. In this study, we changed the aspect ratio of a particle by lengthening the particle and kept its diameter unchanged. Therefore, the E e k s of longer particles (with larger masses) are smaller.
air along parabolic trajectories. The ejected particles were moving upward along the parabolic trajectory before 20 ms, and their velocities were hence decreasing, resulting in the decreasing kinetic energy. It is easy to find that the peak values of s for spherical particles are much larger than those of elongated particles for both glass and brass impactors because spherical particles roll more easily than elongated ones. Comparatively, the s of different elongated particles lengths are obviously smaller than those of spherical particles, and s are slightly smaller as the particles become longer, whether for the glass or brass impactor. We can suppose that the kinetic energy transferred to the ejected particles is to support their movement and more energy will be needed to drive a particle moving if its mass becomes larger. In this study, we changed the aspect ratio of a particle by lengthening the particle and kept its diameter unchanged. Therefore, the s of longer particles (with larger masses) are smaller.

The Impactor
The responses of the impactor are investigated by DEM simulations here, including the penetration depth of the impactor (H0) and the energy of the impactor ( ). Figure 6 shows the normalized penetration depth of granular beds for spheres and elongated particles. Comparing the H0/D0 of brass and glass impactors in Figure 6a,b, it can be found that the final H0 is negatively correlated to the particle length of the granular bed. For the brass impactor, the largest H0 is about 1.5 D0 for the granular bed made of spheres, and the granular beds with Lp = 10 mm elongated particles have the smallest H0 of about 1.1 D0. The shape of the H0(t)/D0 function of the glass impactor is similar to that of the brass impactor, but its value is less than that of the brass impactor. In addition, the glass impactor stopped earlier than the brass impactor, indicating that the brass impactor with larger kinetic energy took more time to halt the penetration. The time of static state of the glass impactor is earlier than that of the brass impactor, and the static state is defined as the unchanged value of H0 of the impactor.

The Impactor
The responses of the impactor are investigated by DEM simulations here, including the penetration depth of the impactor (H 0 ) and the energy of the impactor (E i k ). Figure 6 shows the normalized penetration depth of granular beds for spheres and elongated particles. Comparing the H 0 /D 0 of brass and glass impactors in Figure 6a,b, it can be found that the final H 0 is negatively correlated to the particle length of the granular bed. For the brass impactor, the largest H 0 is about 1.5 D 0 for the granular bed made of spheres, and the granular beds with L p = 10 mm elongated particles have the smallest H 0 of about 1.1 D 0 . The shape of the H 0 (t)/D 0 function of the glass impactor is similar to that of the brass impactor, but its value is less than that of the brass impactor. In addition, the glass impactor stopped earlier than the brass impactor, indicating that the brass impactor with larger kinetic energy took more time to halt the penetration. The time of static state of the glass impactor is earlier than that of the brass impactor, and the static state is defined as the unchanged value of H 0 of the impactor. Because of the different penetration depths of the impactor, the initial potential energy is different, resulting in a slight difference of the Ek0 between the glass impactor and brass impactor. As shown in Figure 7a,b, the Ek0 for the glass and brass impactors at t = 0 ms are 98% and 96.5%, respectively, indicating that the initial total energy of the impactors dominated by Ek0 and its Ep0 is very small. It can also be found from Figure 7a,b that the Because of the different penetration depths of the impactor, the initial potential energy is different, resulting in a slight difference of the E k0 between the glass impactor and brass impactor. As shown in Figure 7a,b, the E k0 for the glass and brass impactors at t = 0 ms are 98% and 96.5%, respectively, indicating that the initial total energy of the impactors dominated by E k0 and its E p0 is very small. It can also be found from Figure 7a,b that the kinetic energy of the impactor E i k decreases rapidly from 0 ms to 5 ms. For both glass and brass impactors, the decrease of E i k slightly increases with the increasing particle length, indicating that the dissipation of E i k is faster with a larger particle length. Furthermore, we quantify the average E i k of five lengths, and it shows that the average E i k obeys an exponential-like dissipation. In addition, the brass impactor dissipates 95% of E i k in more than 10 ms, while the glass impactor dissipates 95% of E i k in less than 5 ms. The main reason is that the initial total energy of the glass impactor is smaller than that of the brass impactor. Because of the different penetration depths of the impactor, the initial potential energy is different, resulting in a slight difference of the Ek0 between the glass impactor and brass impactor. As shown in Figure 7a,b, the Ek0 for the glass and brass impactors at t = 0 ms are 98% and 96.5%, respectively, indicating that the initial total energy of the impactors dominated by Ek0 and its Ep0 is very small. It can also be found from Figure 7a,b that the kinetic energy of the impactor decreases rapidly from 0 ms to 5 ms. For both glass and brass impactors, the decrease of slightly increases with the increasing particle length, indicating that the dissipation of is faster with a larger particle length. Furthermore, we quantify the average of five lengths, and it shows that the average obeys an exponential-like dissipation. In addition, the brass impactor dissipates 95% of in more than 10 ms, while the glass impactor dissipates 95% of in less than 5 ms. The main reason is that the initial total energy of the glass impactor is smaller than that of the brass impactor.

The Particles in the Cup
To study the dynamic response of the particles in the cup, the average contact force of the granular bed (F c pp ), the solid volume fraction (φ p ), the granular temperature (T p ), and the kinetic energy are discussed.
The correlation between H 0 and L p can be explained by the average contact force between particles F c pp in the granular bed. F c pp is the average resultant force of the normal and tangential forces at all contact points, which can describe the strength to resist the penetration of the impactor, as shown in Figure 8. When the impactor penetrates the granular bed, the particles in the cup move randomly and intensively in the first 10 ms. During this period, the impactor and the particles in the granular bed will be in contact and separated from each other until most of the energy of the impactor is dissipated. In Figure 7, it is easy to find that the energy of the impactor sharply decreases in the first 10 ms. The contact forces show a zig-zagging change, which is a characteristic of the discontinuous medium after impact. The contact force is transmitted through the force chain. The particles are discrete, and the force chain structure in the granular system will change at every moment due to the frequent formation and breakage of the contact forces between particles, which leads to the zig-zagging change of the contact force. Furthermore, a longer particle can easily trigger stronger contact force between particles and has a stronger resistance to the penetration of the impactor, thus reducing the penetration depth of the impactor. This can be used to explain the phenomenon shown in Figure 6. cles are discrete, and the force chain structure in the granular system will change at every moment due to the frequent formation and breakage of the contact forces between particles, which leads to the zig-zagging change of the contact force. Furthermore, a longer particle can easily trigger stronger contact force between particles and has a stronger resistance to the penetration of the impactor, thus reducing the penetration depth of the impactor. This can be used to explain the phenomenon shown in Figure 6. For the analysis of granular temperature, solid volume fraction, and kinetic energy, the granular bed was divided into three cylindrical regions to investigate the dynamic responses of the particles in the cup, as shown in Figure 9. To distinguish the particles right under the impactor and surrounding the impactor, the radius of region Ⅰ is chosen to be close to the radius of the impactor. Meanwhile, the thickness of region Ⅲ cannot be smaller than the length of the longest particle (~2/7Rc). Therefore, the radial widths of region Ⅰ, Ⅱ, and Ⅲ are 2/7Rc, 3/7Rc, and 2/7Rc, respectively. The edge between two adjacent regions is fixed, and particles are free to enter or leave one region. For the analysis of granular temperature, solid volume fraction, and kinetic energy, the granular bed was divided into three cylindrical regions to investigate the dynamic responses of the particles in the cup, as shown in Figure 9. To distinguish the particles right under the impactor and surrounding the impactor, the radius of region I is chosen to be close to the radius of the impactor. Meanwhile, the thickness of region III cannot be smaller than the length of the longest particle (~2/7R c ). Therefore, the radial widths of region I, II, and III are 2/7R c , 3/7R c , and 2/7R c , respectively. The edge between two adjacent regions is fixed, and particles are free to enter or leave one region.    Figure 10b. It can be found that the ϕps of region Ⅰ and Ⅱ are closer, and those of region Ⅲ are slightly smaller. It is interesting to find that the ϕp of the granular bed and the H0 of impactor both decrease with longer particles, which is different from previous studies on spherical particles [13,17]. Due to the random distribution of elongated particles, a 3D cage structure will be formed, and there will be more voids in the structure compared with the granular bed of spheres. Therefore, a smaller ϕp of the bed will be formed by the random packing of longer elongated particles. Figure 9. Diagram of sub-region division for granular bed. Region I: 0 < R 1 < 2/7R c ; region II: 2/7R c < R 2 < 5/7R c ; and region III: 5/7R c < R 3 < R c . (The particles count in one region if their mass centers fall in this region; R c is the radius of the cup). Figure 10a shows the comparison of solid volume fractions for granular beds composed of spheres, L p = 4 mm, 6 mm elongated particles, impacted by glass and brass impactors in simulations and experiments. The initial solid volume fractions φ p s of the granular beds are compared in Figure 10b. It can be found that the φ p s of region I and II are closer, and those of region III are slightly smaller. It is interesting to find that the φ p of the granular bed and the H 0 of impactor both decrease with longer particles, which is different from previous studies on spherical particles [13,17]. Due to the random distribution of elongated particles, a 3D cage structure will be formed, and there will be more voids in the structure compared with the granular bed of spheres. Therefore, a smaller φ p of the bed will be formed by the random packing of longer elongated particles. Figure 10a shows the comparison of solid volume fractions for granular beds composed of spheres, Lp = 4 mm, 6 mm elongated particles, impacted by glass and brass impactors in simulations and experiments. The initial solid volume fractions ϕps of the granular beds are compared in Figure 10b. It can be found that the ϕps of region Ⅰ and Ⅱ are closer, and those of region Ⅲ are slightly smaller. It is interesting to find that the ϕp of the granular bed and the H0 of impactor both decrease with longer particles, which is different from previous studies on spherical particles [13,17]. Due to the random distribution of elongated particles, a 3D cage structure will be formed, and there will be more voids in the structure compared with the granular bed of spheres. Therefore, a smaller ϕp of the bed will be formed by the random packing of longer elongated particles.
(a) (b) Figure 10. (a) Comparison of solid volume fraction ϕp as a function of particle length Lp between experiments and simulations for both glass and brass impactors. (b) Initial solid volume fraction ϕp of the whole granular bed or three regions of granular bed as a function of particle length Lp before impact. The particle length of 2 mm represents the data of 2 mm spheres here.
Granular temperature [38] Tp is calculated to describe the fluctuation of translational velocities of the elongated particles in the cup, as shown in Figure 11. The granular temperature of all particles in the cup can be defined as follows: Figure 10. (a) Comparison of solid volume fraction φ p as a function of particle length L p between experiments and simulations for both glass and brass impactors. (b) Initial solid volume fraction φ p of the whole granular bed or three regions of granular bed as a function of particle length L p before impact. The particle length of 2 mm represents the data of 2 mm spheres here.
Granular temperature [38] T p is calculated to describe the fluctuation of translational velocities of the elongated particles in the cup, as shown in Figure 11. The granular temperature of all particles in the cup can be defined as follows: T p = (T p,x +T p,y +T p,z )/3 (11) in which T p,x = (u x − u x ) 2 , T p,y = u y − u y 2 , T p,z = (u z − u z ) 2 , and u − u is the fluctuating velocities of the particles. The operator u is used to calculate the average velocity over all particles. The T p s of three regions are calculated respectively. = (T p,x + T p,y + T p,z ) / 3 (11) in which T p,x = 〈( − 〈 〉)〉 , T p,y = 〈 − 〈 〉 〉 , T p,z = 〈( − 〈 〉)〉 , and − 〈 〉 is the fluctuating velocities of the particles. The operator 〈 〉 is used to calculate the average velocity over all particles. The Tps of three regions are calculated respectively. In Figure 11, it can be found that the difference between Tps and Lps is limited, indicating that Lp has no significant effect on particle movement. The particle temperature of region Ⅰ experiences its peak value in the first 1 ms, and the other two regions are about 3 ms and 6 ms, respectively, which shows that the particles move from the middle to the periphery. The peak value of Tp of the brass impactor is about twice that of the glass impactor in region Ⅰ and region Ⅱ. In short, the granular temperature in region Ⅰ experiences the most dramatic change, followed by region Ⅱ and region Ⅲ.  Figure 12a,b shows the kinetic energy of particles in the cup , which is similar to granular temperature. The particle length has little effect on the dissipation of for any region. The particles left in the cup did not move violently but slightly; therefore, the kinetic energy of the particles left in the cup will not be affected by the particle length. In addition, particle length has little effect on the dissipation of for any region. Comparing the energy dissipation of the three regions, it can be found that the of region Ⅰ has the most drastic change, while that of region Ⅲ has the smallest change. In addition, the In Figure 11, it can be found that the difference between T p s and L p s is limited, indicating that L p has no significant effect on particle movement. The particle temperature of region I experiences its peak value in the first 1 ms, and the other two regions are about 3 ms and 6 ms, respectively, which shows that the particles move from the middle to the periphery. The peak value of T p of the brass impactor is about twice that of the glass impactor in region I and region II. In short, the granular temperature in region I experiences the most dramatic change, followed by region II and region III. Figure 12a,b shows the kinetic energy of particles in the cup E c k , which is similar to granular temperature. The particle length has little effect on the dissipation of E c k for any region. The particles left in the cup did not move violently but slightly; therefore, the kinetic energy of the particles left in the cup will not be affected by the particle length. In addition, particle length has little effect on the dissipation of E c k for any region. Comparing the energy dissipation of the three regions, it can be found that the E c k of region I has the most drastic change, while that of region III has the smallest change. In addition, the total E c k of three regions of glass impactor is larger than that of brass impactor, which indicates that the impactor with less initial energy transfers its higher energy to the particles in the cup. The brass impactor with higher initial energy transfers more energy to the ejected particles, as shown in Figure 5b.

Effect of Friction Coefficient
To evaluate the influence of particle friction on penetration, the friction coefficients of the target particles (μp-p), impactor and particles (μi-p), and cup wall and particles (μw-p) were studied. We chose the steel elongated particles of Lp = 6 mm with initial friction coefficients and the brass impactor for the DEM simulation in this section. Figure 13 shows that the effect of μp-p is dominant. The ejecta mass decreases dramatically when μp-p ranges from 0 to 0.2, and it decreases gently as μp-p changes from 0.2 to 1.0 gradually, while μi-p and μw-p have very little effect on ejecta mass. Figure 13. Relationship between percentage of ejecta mass and friction coefficient μp-p, μi-p, and μw-p

Effect of Friction Coefficient
To evaluate the influence of particle friction on penetration, the friction coefficients of the target particles (µ p-p ), impactor and particles (µ i-p ), and cup wall and particles (µ w-p ) were studied. We chose the steel elongated particles of L p = 6 mm with initial friction coefficients and the brass impactor for the DEM simulation in this section. Figure 13 shows that the effect of µ p-p is dominant. The ejecta mass decreases dramatically when µ p-p ranges from 0 to 0.2, and it decreases gently as µ p-p changes from 0.2 to 1.0 gradually, while µ i-p and µ w-p have very little effect on ejecta mass.
The non-dimensional penetration depth (H 0 /D 0 ) is used to illustrate the correlation among various friction coefficients between particles, which is exhibited in detail in Figure 14a. It is interesting to find that the maximum H 0 of the impactor is about 2D 0 when µ p-p = 0 for condition I, and H 0 decreases with increasing µ p-p. The impactor will rebound to the upside and not penetrate into the granular bed while µ p-p ≥ 0.7. However, the friction coefficients of µ i-p and µ w-p have little effect on H 0 . The maximum values of H 0 for different µ i-p s are less than 0.5 D 0 , as shown in Figure 14b, and the H 0 hardly changes for various µ w-p s in Figure 14c. The vertical velocity of impactor V is shown in Figure 15a,b. When µ p-p is less than 0.7, the decrease of V is positively correlated to µ p-p , as shown in Figure 15a. The changes of V for different µ i-p s and µ w-p s are similar, as shown in Figure 15b.
of the target particles (μp-p), impactor and particles (μi-p), and cup wall and particles (μw-p) were studied. We chose the steel elongated particles of Lp = 6 mm with initial friction coefficients and the brass impactor for the DEM simulation in this section. Figure 13 shows that the effect of μp-p is dominant. The ejecta mass decreases dramatically when μp-p ranges from 0 to 0.2, and it decreases gently as μp-p changes from 0.2 to 1.0 gradually, while μi-p and μw-p have very little effect on ejecta mass. The non-dimensional penetration depth (H0/D0) is used to illustrate the correlation among various friction coefficients between particles, which is exhibited in detail in Figure  14a. It is interesting to find that the maximum H0 of the impactor is about 2D0 when μp-p = 0 for condition Ⅰ, and H0 decreases with increasing μp-p. The impactor will rebound to the upside and not penetrate into the granular bed while μp-p ≥ 0.7. However, the friction coefficients of μi-p and μw-p have little effect on H0. The maximum values of H0 for different μi-ps are less than 0.5 D0, as shown in Figure 14b, and the H0 hardly changes for various μw-ps in Figure 14c. The vertical velocity of impactor V is shown in Figure 15a,b. When μpp is less than 0.7, the decrease of V is positively correlated to μp-p, as shown in Figure 15a. The changes of V for different μi-ps and μw-ps are similar, as shown in Figure 15b.  Figure 16 shows the velocity profiles of particles with μp-p = 0.8 at different times im pacted by the brass impactor. The color of the granular bed represents the velocity ma nitude, which is mainly in the range of 0.3 m/s-1 m/s. From Figures 14a and 16, it can b found that the impactor rebounds to the height of 0.5D0 at t = 50 ms. Then, the impact starts falling under gravity, and its velocity reaches a small peak at t = 90 ms, as shown Figure 15a. Some ejected particles will contact with the impactor during t = 90-100 ms, shown in Figure 16, which slows the impactor slightly. The second penetration begins t = 90 ms, corresponding to the decreasing velocity, and the contact position is higher b cause the surface of the granular bed is irregular after the first impact. The maximum H is kept at 0.15D0 when t = 150 ms, indicating the termination of the whole penetration. A a whole, the velocity direction of the impactor is not always downward during its pen tration, resulting in the spatial curve of the actual penetration path of the impactor. addition, when the impactor rebounds, the trajectory of the impactor is a parabola, whic makes the drop-point of the second penetration depart from that of the first one.  Figure 16 shows the velocity profiles of particles with µ p-p = 0.8 at different times impacted by the brass impactor. The color of the granular bed represents the velocity magnitude, which is mainly in the range of 0.3 m/s-1 m/s. From Figures 14a and 16, it can be found that the impactor rebounds to the height of 0.5D 0 at t = 50 ms. Then, the impactor starts falling under gravity, and its velocity reaches a small peak at t = 90 ms, as shown in Figure 15a. Some ejected particles will contact with the impactor during t = 90-100 ms, as shown in Figure 16, which slows the impactor slightly. The second penetration begins at t = 90 ms, corresponding to the decreasing velocity, and the contact position is higher because the surface of the granular bed is irregular after the first impact. The maximum H 0 is kept at 0.15D 0 when t = 150 ms, indicating the termination of the whole penetration. As a whole, the velocity direction of the impactor is not always downward during its penetration, resulting in the spatial curve of the actual penetration path of the impactor. In addition, when the impactor rebounds, the trajectory of the impactor is a parabola, which makes the drop-point of the second penetration depart from that of the first one. Comparing the ejecta mass and penetration depth under different particle frictions (Figures 13 and 14a), we presume that the objective particles would contact each other with a large contact force and build a strong quasi-continuous medium under a strong force chain to resist the penetration of the impactor, resulting in the decrease of V, as Comparing the ejecta mass and penetration depth under different particle frictions (Figures 13 and 14a), we presume that the objective particles would contact each other with a large contact force and build a strong quasi-continuous medium under a strong force chain to resist the penetration of the impactor, resulting in the decrease of V, as shown in Figure 15a. Therefore, we can say that the ejecta mass for a larger friction coefficient is decreased by the stronger force chain, because it takes more energy to separate the particles. Based on this consideration, we study the vertical contact force between the impactor and particles F c ip , and the average contact force between particles and particles F c pp , as shown in Figure 17a,b, respectively. It is obvious that the forces F c ip and F c pp increase with increasing µ p-p , and the force chain is the strongest when µ p-p = 1. Figure 16. Velocity profiles of elongated particles (Lp = 6 mm) impacted by brass impactor. (Lp = 6 mm, μp-p = 0.8, arrow represents velocity of impactor, whose color remains unchanged.) Comparing the ejecta mass and penetration depth under different particle frictions (Figures 13 and 14a), we presume that the objective particles would contact each other with a large contact force and build a strong quasi-continuous medium under a strong force chain to resist the penetration of the impactor, resulting in the decrease of V, as shown in Figure 15a. Therefore, we can say that the ejecta mass for a larger friction coefficient is decreased by the stronger force chain, because it takes more energy to separate the particles. Based on this consideration, we study the vertical contact force between the impactor and particles , and the average contact force between particles and particles , as shown in Figure 17a,b, respectively. It is obvious that the forces and increase with increasing μp-p, and the force chain is the strongest when μp-p = 1.

Effect of Particle Configuration
To find the effects of the shape of the cup and particle configuration, we chose steel elongated particles of L p = 10 mm and the brass impactor to pack in the cylindrical and cuboid cup in the lateral and vertical directions, respectively, as shown in Figure 18. The friction coefficients and material parameters are shown in Table 2. In order to make the elongated particles more regular when they are arranged horizontally, the size of the cuboid cup is the integral multiple of the length of the elongated particles. The solid volume fractions φ p s are shown in Figure 19. It can be seen that the φ p s of the cuboid cup are larger than that of the cylindrical cup. The φ p of vertical arrangement is the largest, and the φ p of random packing is the smallest.
The orientational order of particles can be monitored by diagonalization of the symmetric traceless order tensor Q [38], as follows: where l n i(j) is the unit vector along the major axis of the elongated particle n among all N particles. The largest eigenvalue of Q is the primary order parameter S r , which quantifies the degree of alignment. When i = j, δ ij = 1, or δ ij = 0, the order parameter S r is equal to one if all particles are packed in the same direction, and S r is zero if the orientation of every particle is different. friction coefficients and material parameters are shown in Table 2. In order to make the elongated particles more regular when they are arranged horizontally, the size of the cuboid cup is the integral multiple of the length of the elongated particles. The solid volume fractions ϕps are shown in Figure 19. It can be seen that the ϕps of the cuboid cup are larger than that of the cylindrical cup. The ϕp of vertical arrangement is the largest, and the ϕp of random packing is the smallest.  The orientational order of particles can be monitored by diagonalization of the symmetric traceless order tensor Q [38], as follows: (12) Figure 19. Comparison between solid volume fractions for different packing types (type 1-lateral arrangement, type 2-vertical arrangement, and type 3-random arrangement). The length of particles in the cup is 10 mm.
Although S r can describe the order of particle configuration, it cannot show the direction of particle orientation. For example, if the major axes of all particles are all in the same inclined direction, the direction of particles cannot be specifically distinguished using Equation (12). Therefore, the average orientational parameter O of all particles is used as the quantitative analysis of the particle orientation. The orientational parameter [39] of a single elongated particle is calculated by the acute angle between the major axis of the elongated particle and the vertical axis. A larger O leads to a more vertical particle. Figure 20a,b show the order parameters S r and orientational parameters O for three packing types in the cylindrical and cuboid cups. The values of S r in the cuboid cup are all bigger than those in the cylindrical cup, which is caused by the restriction of the cylindrical wall.
The orientational order of particles can be monitored by diagonalization of the sym-metric traceless order tensor Q [38], as follows: where ( ) is the unit vector along the major axis of the elongated particle among all N particles. The largest eigenvalue of is the primary order parameter Sr, which quantifies the degree of alignment. When i = j, δ ij = 1, or δ ij = 0, the order parameter Sr is equal to one if all particles are packed in the same direction, and Sr is zero if the orientation of every particle is different.
Although Sr can describe the order of particle configuration, it cannot show the direction of particle orientation. For example, if the major axes of all particles are all in the same inclined direction, the direction of particles cannot be specifically distinguished using Equation (12). Therefore, the average orientational parameter O of all particles is used as the quantitative analysis of the particle orientation. The orientational parameter [39] of a single elongated particle is calculated by the acute angle between the major axis of the elongated particle and the vertical axis. A larger O leads to a more vertical particle. Figure  20a,b show the order parameters Sr and orientational parameters O for three packing types in the cylindrical and cuboid cups. The values of Sr in the cuboid cup are all bigger than those in the cylindrical cup, which is caused by the restriction of the cylindrical wall.
(a) (b) Figure 20. The order parameters and orientational parameters of three packings in cylindrical cup and cuboid cup: (a) is the order parameters for three packing types; (b) is the orientational parameters for three packing types (type 1-lateral arrangement, type 2-vertical arrangement, and type 3-random arrangement). The length of particles in the cup is L p = 10 mm.
As shown in Figure 21, the ejecta masses of the cylindrical cup are all larger than that of the cuboid cup. The ejecta mass of lateral arrangement is the largest, and that of vertical arrangement is the smallest. The elongated particles of lateral arrangement will easily roll out of the cup from the cup edge during impact, which leads to the larger ejecta masses of lateral arrangement. As shown in Figure 22, the H 0 of the random packing is the largest, and the H 0 of the vertical packing is the smallest in the two cups. In short, the cuboid cup has larger S r and φ p , which lead to the smaller ejecta mass and penetration depth of the cuboid cup compared to those of the cylindrical cup. In addition, the vertical arrangement of elongated particles is more regular and has a stronger ability to resist penetration than the lateral arrangement. As shown in Figure 21, the ejecta masses of the cylindrical cup are all larger than that of the cuboid cup. The ejecta mass of lateral arrangement is the largest, and that of vertical arrangement is the smallest. The elongated particles of lateral arrangement will easily roll out of the cup from the cup edge during impact, which leads to the larger ejecta masses of lateral arrangement. As shown in Figure 22, the H0 of the random packing is the largest, and the H0 of the vertical packing is the smallest in the two cups. In short, the cuboid cup has larger Sr and ϕp, which lead to the smaller ejecta mass and penetration depth of the cuboid cup compared to those of the cylindrical cup. In addition, the vertical arrangement of elongated particles is more regular and has a stronger ability to resist penetration than the lateral arrangement.

Conclusions
Three-dimensional DEM simulations and experiments of spheres impacting into elongated particles at a low speed are investigated in this study. By validating the numerical simulations with the experimental results, we demonstrate the feasibility of using the

Conclusions
Three-dimensional DEM simulations and experiments of spheres impacting into elongated particles at a low speed are investigated in this study. By validating the numerical simulations with the experimental results, we demonstrate the feasibility of using the above code to understand the sphere impacting elongated particles. Close agreement between the DEM simulations and experimental results can be obtained in terms of ejecta masses and initial solid volume fraction. The effects of the particle length, the friction, and the configuration of the elongated particle bed on impact are discussed in this study.
The effect of the sphere and elongated particles of four particles lengths on penetration depth, ejecta mass, granular temperature, and contact force are studied. Then, the friction coefficients µ p-p , µ i-p , and µ w-p ranging from 0 to 1 are discussed; the impactor rebounds when µ p-p is larger than 0.7. As for the configuration of the elongated particle bed, three packing types and two different cups are compared, ejecta mass of particles and penetration depth of the impactor are studied. Based on the present studies, the following conclusions can be drawn: 1.
The effect of particle length. The ejecta mass of the spherical particle bed is obviously larger than that of the elongated particle bed. The granular bed of longer particles has a smaller penetration depth due to the spatial structure of elongated particles, although the solid volume fraction is smaller. In addition, the average contact force between particles is positively correlated to particle length. The average kinetic energy of the impactor obeys an exponential-like dissipation, and the particle length of the elongated particles has little effect on the energy allocation from the impactor to the ejected particles and particles in the cup.

2.
The effect of friction. The µ p-p has a significant effect on the ejecta mass and penetration depth of the impactor, while µ i-p and µ w-p have a limited effect. The ejecta mass and penetration depth are negatively correlated to µ p-p . The contact force between particles and particles or impactors are positively correlated to µ p-p . 3.
The effect of particle configuration. The cuboid cup can obtain a more dense and regular granular bed. The ejecta mass and penetration depth of vertical arrangement are the smallest. For the same arrangement of elongated particles, the penetration depth is negatively correlated to order parameters and solid volume fraction.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
F i external forces of particle i (N) T i external torques of particle i (N·m) m i mass of particle i (g) v i translational velocity of particle i (m/s) I i moment of inertia of particle i (kg·m 2 ) ω i rotational velocity of the particle i (rad/s) δ n overlap in normal direction (mm) F n normal contact force (N) E * equivalent Young's modulus R * equivalent radius of two objects φ p solid volume fraction F c pp average contact force between particles (N) µ p-p coefficient of friction between target particles µ w-p coefficient of friction between wall and particles µ i-p coefficient of friction between impactor and particles T p,q granular temperature of particles along q (q = x, y, z) direction (m 2 /s 2 ) F c ip vertical contact force between impactor and particles (N) O orientational parameter S r the order parameter